{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 平均脸\n",
    "![这里写图片描述](http://www.learnopencv.com/wp-content/uploads/2016/05/average_entrepreneur.jpg)\n",
    "![这里写图片描述](https://www.learnopencv.com/wp-content/uploads/2016/05/average_best_actress-300x300.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "py代码以及相关数据地址：https://www.learnopencv.com/wp-content/uploads/2016/05/FaceAverage.zip\n",
    "\n",
    "最初博文地址：https://www.learnopencv.com/average-face-opencv-c-python-tutorial/\n",
    "中文翻译：http://blog.csdn.net/GraceDD/article/details/51382952\n",
    "\n",
    "\n",
    "中文改编地址：[《手把手：用OpenCV亲手给小扎、Musk等科技大佬们做一张“平均脸”（附Python代码）》](https://mp.weixin.qq.com/s?__biz=MjM5MTQzNzU2NA==&mid=2651654758&idx=1&sn=b60e2da0b4e9cffed660f44bd624eb9e&chksm=bd4c2df58a3ba4e32f938df33cdc780bd7041087c6f3b82cf059c036a50c12e97067a8d12815&mpshare=1&scene=1&srcid=1123eFDjNTtDFdq4GS8M2e8d#rd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 平均脸\n",
    "在完成各个library的安装后。\n",
    "第一步：将要平均的照片放入faces文档，确保图片为jpg格式。\n",
    "第二步：在终端运行 python face_landmark_detection.py shape_predictor_68_face_landmarks.dat faces，并在程序运行结束后将所有faces文档中的文件复制到presidents文档中（如无法完成dlib安装，可略过该步骤，直接用文摘菌提供的素材）\n",
    "第三步：在终端运行 python faceAverage.py \n",
    "这样就能看到制作成功的平均脸了！\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import cv2\n",
    "import numpy as np\n",
    "import math\n",
    "import sys\n",
    "\n",
    "# Read points from text files in directory\n",
    "def readPoints(path) :\n",
    "    # Create an array of array of points.\n",
    "    pointsArray = [];\n",
    "\n",
    "    #List all files in the directory and read points from text files one by one\n",
    "    for filePath in os.listdir(path):\n",
    "        \n",
    "        if filePath.endswith(\".txt\"):\n",
    "            \n",
    "            #Create an array of points.\n",
    "            points = [];            \n",
    "            \n",
    "            # Read points from filePath\n",
    "            with open(os.path.join(path, filePath)) as file :\n",
    "                for line in file :\n",
    "                    x, y = line.split()\n",
    "                    points.append((int(x), int(y)))\n",
    "            \n",
    "            # Store array of points\n",
    "            pointsArray.append(points)\n",
    "            \n",
    "    return pointsArray;\n",
    "\n",
    "# Read all jpg images in folder.\n",
    "def readImages(path) :\n",
    "    \n",
    "    #Create array of array of images.\n",
    "    imagesArray = [];\n",
    "    \n",
    "    #List all files in the directory and read points from text files one by one\n",
    "    for filePath in os.listdir(path):\n",
    "        if filePath.endswith(\".jpg\"):\n",
    "            # Read image found.\n",
    "            img = cv2.imread(os.path.join(path,filePath));\n",
    "\n",
    "            # Convert to floating point\n",
    "            img = np.float32(img)/255.0;\n",
    "\n",
    "            # Add to array of images\n",
    "            imagesArray.append(img);\n",
    "            \n",
    "    return imagesArray;\n",
    "                \n",
    "# Compute similarity transform given two sets of two points.\n",
    "# OpenCV requires 3 pairs of corresponding points.\n",
    "# We are faking the third one.\n",
    "\n",
    "def similarityTransform(inPoints, outPoints) :\n",
    "    s60 = math.sin(60*math.pi/180);\n",
    "    c60 = math.cos(60*math.pi/180);  \n",
    "  \n",
    "    inPts = np.copy(inPoints).tolist();\n",
    "    outPts = np.copy(outPoints).tolist();\n",
    "    \n",
    "    xin = c60*(inPts[0][0] - inPts[1][0]) - s60*(inPts[0][1] - inPts[1][1]) + inPts[1][0];\n",
    "    yin = s60*(inPts[0][0] - inPts[1][0]) + c60*(inPts[0][1] - inPts[1][1]) + inPts[1][1];\n",
    "    \n",
    "    inPts.append([np.int(xin), np.int(yin)]);\n",
    "    \n",
    "    xout = c60*(outPts[0][0] - outPts[1][0]) - s60*(outPts[0][1] - outPts[1][1]) + outPts[1][0];\n",
    "    yout = s60*(outPts[0][0] - outPts[1][0]) + c60*(outPts[0][1] - outPts[1][1]) + outPts[1][1];\n",
    "    \n",
    "    outPts.append([np.int(xout), np.int(yout)]);\n",
    "    \n",
    "    tform = cv2.estimateRigidTransform(np.array([inPts]), np.array([outPts]), False);  # 多个二维点对之间的仿射变换矩阵（使用误差最小准则）\n",
    "    \n",
    "    return tform;\n",
    "\n",
    "\n",
    "# Check if a point is inside a rectangle\n",
    "def rectContains(rect, point) :\n",
    "    if point[0] < rect[0] :\n",
    "        return False\n",
    "    elif point[1] < rect[1] :\n",
    "        return False\n",
    "    elif point[0] > rect[2] :\n",
    "        return False\n",
    "    elif point[1] > rect[3] :\n",
    "        return False\n",
    "    return True\n",
    "\n",
    "# Calculate delanauy triangle\n",
    "def calculateDelaunayTriangles(rect, points):\n",
    "    # Create subdiv\n",
    "    subdiv = cv2.Subdiv2D(rect);\n",
    "   \n",
    "    # Insert points into subdiv\n",
    "    for p in points:\n",
    "        subdiv.insert((p[0], p[1]));\n",
    "\n",
    "   \n",
    "    # List of triangles. Each triangle is a list of 3 points ( 6 numbers )\n",
    "    triangleList = subdiv.getTriangleList();\n",
    "\n",
    "    # Find the indices of triangles in the points array\n",
    "\n",
    "    delaunayTri = []\n",
    "    \n",
    "    for t in triangleList:\n",
    "        pt = []\n",
    "        pt.append((t[0], t[1]))\n",
    "        pt.append((t[2], t[3]))\n",
    "        pt.append((t[4], t[5]))\n",
    "        \n",
    "        pt1 = (t[0], t[1])\n",
    "        pt2 = (t[2], t[3])\n",
    "        pt3 = (t[4], t[5])        \n",
    "        \n",
    "        if rectContains(rect, pt1) and rectContains(rect, pt2) and rectContains(rect, pt3):\n",
    "            ind = []\n",
    "            for j in range(0, 3):\n",
    "                for k in range(0, len(points)):                    \n",
    "                    if(abs(pt[j][0] - points[k][0]) < 1.0 and abs(pt[j][1] - points[k][1]) < 1.0):\n",
    "                        ind.append(k)                            \n",
    "            if len(ind) == 3:                                                \n",
    "                delaunayTri.append((ind[0], ind[1], ind[2]))\n",
    "\n",
    "    return delaunayTri\n",
    "\n",
    "\n",
    "def constrainPoint(p, w, h) :\n",
    "    p =  ( min( max( p[0], 0 ) , w - 1 ) , min( max( p[1], 0 ) , h - 1 ) )\n",
    "    return p;\n",
    "\n",
    "# Apply affine transform calculated using srcTri and dstTri to src and\n",
    "# output an image of size.\n",
    "def applyAffineTransform(src, srcTri, dstTri, size) :\n",
    "    \n",
    "    # Given a pair of triangles, find the affine transform.\n",
    "    warpMat = cv2.getAffineTransform( np.float32(srcTri), np.float32(dstTri) )\n",
    "    \n",
    "    # Apply the Affine Transform just found to the src image\n",
    "    dst = cv2.warpAffine( src, warpMat, (size[0], size[1]), None, flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REFLECT_101 )\n",
    "\n",
    "    return dst\n",
    "\n",
    "\n",
    "# Warps and alpha blends triangular regions from img1 and img2 to img\n",
    "def warpTriangle(img1, img2, t1, t2) :\n",
    "\n",
    "    # Find bounding rectangle for each triangle\n",
    "    r1 = cv2.boundingRect(np.float32([t1]))\n",
    "    r2 = cv2.boundingRect(np.float32([t2]))\n",
    "\n",
    "    # Offset points by left top corner of the respective rectangles\n",
    "    t1Rect = [] \n",
    "    t2Rect = []\n",
    "    t2RectInt = []\n",
    "\n",
    "    for i in range(0, 3):\n",
    "        t1Rect.append(((t1[i][0] - r1[0]),(t1[i][1] - r1[1])))\n",
    "        t2Rect.append(((t2[i][0] - r2[0]),(t2[i][1] - r2[1])))\n",
    "        t2RectInt.append(((t2[i][0] - r2[0]),(t2[i][1] - r2[1])))\n",
    "\n",
    "\n",
    "    # Get mask by filling triangle\n",
    "    mask = np.zeros((r2[3], r2[2], 3), dtype = np.float32)\n",
    "    cv2.fillConvexPoly(mask, np.int32(t2RectInt), (1.0, 1.0, 1.0), 16, 0);\n",
    "\n",
    "    # Apply warpImage to small rectangular patches\n",
    "    img1Rect = img1[r1[1]:r1[1] + r1[3], r1[0]:r1[0] + r1[2]]\n",
    "    \n",
    "    size = (r2[2], r2[3])\n",
    "\n",
    "    img2Rect = applyAffineTransform(img1Rect, t1Rect, t2Rect, size)\n",
    "    \n",
    "    img2Rect = img2Rect * mask\n",
    "\n",
    "    # Copy triangular region of the rectangular patch to the output image\n",
    "    img2[r2[1]:r2[1]+r2[3], r2[0]:r2[0]+r2[2]] = img2[r2[1]:r2[1]+r2[3], r2[0]:r2[0]+r2[2]] * ( (1.0, 1.0, 1.0) - mask )\n",
    "     \n",
    "    img2[r2[1]:r2[1]+r2[3], r2[0]:r2[0]+r2[2]] = img2[r2[1]:r2[1]+r2[3], r2[0]:r2[0]+r2[2]] + img2Rect\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "'''\n",
    "1.读入图 + 读入关键点信息  readPoints  readImages \n",
    "2.平均脸的眼角位置（这样其他脸，按照眼睛位置对齐）  eyecornerDst\n",
    "3.新的8个初始边界点 boundaryPts （为了后续做脸谱网络用的）\n",
    "4.设置初始平均脸 pointsAvg （随便找个脸68个关键点 + 8个初始点）\n",
    "5.根据眼睛位置，进行人脸初步对齐\n",
    "6.计算初始平均脸的脸谱网络76点（calculateDelaunayTriangles）\n",
    "7.根据脸谱网络二次人脸对齐\n",
    "\n",
    "'''\n",
    "# 1.读入图 + 读入关键点信息\n",
    "path = 'presidents/'\n",
    "\n",
    "# Dimensions of output image\n",
    "w = 600;\n",
    "h = 600;\n",
    "\n",
    "# Read points for all images\n",
    "allPoints = readPoints(path);\n",
    "allPoints\n",
    "\n",
    "# Read all images\n",
    "images = readImages(path);\n",
    "images\n",
    "\n",
    "# 2.平均脸的眼角位置（这样其他脸，按照眼睛位置对齐）\n",
    "# 对齐的眼角预先设定好\n",
    "# 确保两只眼睛的点都在一个水平线上，面部中心大约在离顶端三分之一高度的位置。所以我将眼角位置设为（0.3*宽，高/3）和（0.7*宽，高/3）。\n",
    "eyecornerDst = [ (np.int(0.3 * w ), np.int(h / 3)), (np.int(0.7 * w ), np.int(h / 3)) ];\n",
    "\n",
    "imagesNorm = [];\n",
    "pointsNorm = [];\n",
    "\n",
    "# 3.新的8个初始边界点 boundaryPts\n",
    "# Add boundary points for delaunay triangulation  根据图像长宽、设定8个初始边界点\n",
    "boundaryPts = np.array([(0,0), (w/2,0), (w-1,0), (w-1,h/2), ( w-1, h-1 ), ( w/2, h-1 ), (0, h-1), (0,h/2) ]);  \n",
    "\n",
    "# 4. 初始平均脸\n",
    "# Initialize location of average points to 0s\n",
    "pointsAvg = np.array([(0,0)]* ( len(allPoints[0]) + len(boundaryPts) ), np.float32());  #　初始化平均脸\n",
    "\n",
    "n = len(allPoints[0]);\n",
    "\n",
    "numImages = len(images)\n",
    "\n",
    "\n",
    "# 5.根据眼睛位置，进行人脸初步对齐\n",
    "# Warp images and trasnform landmarks to output coordinate system,\n",
    "# and find average of transformed landmarks.\n",
    "\n",
    "for i in range(0, numImages):\n",
    "\n",
    "    points1 = allPoints[i];\n",
    "\n",
    "    # Corners of the eye in input image\n",
    "    eyecornerSrc  = [ allPoints[i][36], allPoints[i][45] ] ;\n",
    "\n",
    "    # Compute similarity transform  眼角对齐转换矩阵\n",
    "    tform = similarityTransform(eyecornerSrc, eyecornerDst);　　# 2 * 3 转换矩阵\n",
    "\n",
    "    # Apply similarity transformation　通过眼角对齐转换矩阵 进行人脸对齐\n",
    "    img = cv2.warpAffine(images[i], tform, (w,h));\n",
    "\n",
    "    # Apply similarity transform on points\n",
    "    points2 = np.reshape(np.array(points1), (68,1,2));      # (68, 1, 2)  \n",
    "\n",
    "    points = cv2.transform(points2, tform);   # (68, 1, 2)  * (2,3)  --> (68, 1, 2)\n",
    "\n",
    "    points = np.float32(np.reshape(points, (68, 2)));  # (68, 2)\n",
    "\n",
    "    # Append boundary points. Will be used in Delaunay Triangulation\n",
    "    points = np.append(points, boundaryPts, axis=0)   # point（68个） + boundaryPts (8个边界点)  (76, 2)\n",
    "\n",
    "    # Calculate location of average landmark points.\n",
    "    pointsAvg = pointsAvg + points / numImages;  # (76, 2)  初始平均脸\n",
    "\n",
    "    pointsNorm.append(points);  # 6*72*2\n",
    "    imagesNorm.append(img);  \n",
    "\n",
    "\n",
    "# 6.计算初始平均脸的脸谱网络76点（calculateDelaunayTriangles）\n",
    "# Delaunay triangulation\n",
    "rect = (0, 0, w, h);\n",
    "dt = calculateDelaunayTriangles(rect, np.array(pointsAvg));  # delanauy三角计算 142*3\n",
    "\n",
    "# 7 根据脸谱网络二次人脸对齐\n",
    "# Output image\n",
    "output = np.zeros((h,w,3), np.float32());\n",
    "\n",
    "# Warp input images to average image landmarks\n",
    "for i in range(0, len(imagesNorm)) :\n",
    "    img = np.zeros((h,w,3), np.float32());\n",
    "    # Transform triangles one by one\n",
    "    for j in range(0, len(dt)) :\n",
    "        tin = []; \n",
    "        tout = [];\n",
    "\n",
    "        for k in range(0, 3) :                \n",
    "            pIn = pointsNorm[i][dt[j][k]];     #需要对齐的脸的脸谱网络计算\n",
    "            pIn = constrainPoint(pIn, w, h);  #　约束wh，如果有越界的，要归到图像内\n",
    "\n",
    "            pOut = pointsAvg[dt[j][k]];      # 平均脸的脸谱网络计算\n",
    "            pOut = constrainPoint(pOut, w, h);\n",
    "\n",
    "            tin.append(pIn);  \n",
    "            tout.append(pOut);\n",
    "\n",
    "        warpTriangle(imagesNorm[i], img, tin, tout);  # 把需要对齐的脸谱网络，向平均脸脸谱网络 转换\n",
    "\n",
    "\n",
    "    # Add image intensities for averaging\n",
    "    output = output + img;\n",
    "\n",
    "# Divide by numImages to get average\n",
    "output = output / numImages;\n",
    "\n",
    "# Save result\n",
    "cv2.imwrite('myaverageface.png', (output * 255).astype('uint8'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
